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Abstract.—Phidiana hiltoni is a conspicuous nudibranch sea slug native to the north¬ 
eastern Pacific Ocean. Over the past thirty years the range of P. hiltoni has expanded 
about 200 km northward, but the mechanism that facilitated this expansion is poorly 
understood. In this study, we use mtDNA and microsatellite data to investigate the 
population structure of P hiltoni in its historical range as well as in recently colonized 
localities. Microsatellite analyses reveal little to no genetic structure and thus high gene 
flow throughout the range of P. hiltoni. This is consistent with mtDNA analysis results, 
which revealed shared haplotypes between Southern, Central and Northern popula¬ 
tions. However, AMOVA of mtDNA data did recover some genetic structure among 
geographic regions. This, along with same group memberships in the microsatellite 
data of individuals from sites like Cave Landing, suggest a certain degree of local 
recruitment and reduced vagility. Recently established populations in Northern Cal¬ 
ifornia contain two unique mtDNA haplotypes that are not present elsewhere, but 
microsatellite data do not differentiate these from other populations. The mismatch 
between mtDNA and microsatellite data could be explained by the mating system of 
this aggressive, hermaphroditic species as well as the sporadic nature of the northward 
dispersal. Analyses of historical abundance data of P. hiltoni suggest a population de¬ 
cline in Southern California. Together, these results suggest a northward population 
shift, rather than a range expansion, possibly related to ongoing changes in nearshore 
oceanographic conditions in the region. 


Rising ocean temperatures driven by global climate change are having dramatic impacts 
on coastal ecosystems around the world (McGowan et al. 1998; Sorte et al. 2011). One 
of the most noticeable effects is the poleward range expansion of certain species (Dawson 
et al. 2010; Sorte et al. 2011; Sunday et al. 2012; Canning-Clode and Carlton 2017). Par¬ 
ticularly problematic are range expansions of predatory species, which can have significant 
impacts on the trophic structure of newly colonized ecosystems (Zeidberg and Robinson 
2007; Gallardo et al. 2016). However, not all of these range expansions are permanent; 
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some result from regular oscillations in ocean temperatures (e.g., El Nino events). In these 
cases, populations often return to their original range following ephemeral warming events, 
making it difficult to attribute in the short term any particular range shift to longer term 
climate change (Schultz et al. 2011). Poleward range expansions can also be accompanied 
by extirpation at lower latitudes, resulting in shifts at both ends of species ranges (Parme¬ 
san et al. 1999; Bates et al. 2014). Range shifts may constitute a more pervasive indication 
of permanent changes in the ecological structure of biotas as they can be more difficult to 
reverse (Parmesan et al. 1999; Schultz et al. 2011). However, range shifts are difficult to 
detect and precisely quantify, particularly in marine species with low abundance and/or 
population densities (Bates et al. 2015). 

Phidicma hiltoni is a relatively large and conspicuous aeolid nudibranch native to 
the northeastern Pacific Ocean. Like most nudibranchs, P. hiltoni is a simultaneous 
hermaphrodite, but the mating behavior of this species is poorly understood. The diet 
of P. hiltoni consists mostly of hydroids and other cnidarians; however, individuals of 
this species are known to attack and consume other sea slugs, particularly small, soft- 
bodied aeolids and dendronotaceans, including conspecifics (Goddard et al. 2011). Phidi- 
ana hiltoni has relatively large eggs and lecithotrophic larval development (Goddard 2004); 
its larvae do not need to feed in the water column and are capable of settlement and meta¬ 
morphosis within a day or two of hatching. Thus, compared to planktotrophic species, 
dispersal by the larvae of P. hiltoni is greatly reduced. Historically, P. hiltoni was found as 
far south as Isla Cedros off the coast of Baja California, Mexico and as far north as Pacific 
Grove, California (Goddard et al. 2011). In 1977, Phidicma hiltoni was discovered north of 
Monterey Bay (Goddard et al. 2011). Once across Monterey Bay, P hiltoni rapidly made 
its way up the coast; it was found just north of San Francisco Bay (Duxbury Reef, Marin 
County) in 1992 and now is present as far north as Bodega Bay, California, representing a 
200-km northward range expansion in 40 years (Goddard et al. 2011; Goddard et al. 2018). 
At Duxbury Reef, P. hiltoni quickly became the dominant sea slug, with apparent negative 
impacts on other nudibranch species, likely through a combination of direct predation and 
competition for shared hydroid prey (Goddard et al. 2011). 

The mechanism behind the range expansion of P hiltoni is not well understood, but 
has been potentially linked to warming coastal waters and shifts in ocean currents along 
the California coast (Schultz et al. 2011). Changes in ocean circulation, which drives lar¬ 
val transport, can potentially increase the risk of species introductions and/or dispersals 
(Harley et al. 2006; Sorte et al. 2011; Wilson et al. 2016). Whereas the dispersal potential of 
species with planktonic feeding larvae is relatively well understood (Scheltema 1986), less is 
known about how species with lecithotrophic development may respond to oceanographic 
changes. A majority of lecithotrophic sea slug species are found in warmer, nutrient-poor 
waters, where having a short-lived, non-feeding larval stage can lead to reductions in larval 
mortality at the cost of reduced fecundity and vagility (Goddard 2004; Goddard and Her- 
mosillo 2008). Phidicma hiltoni is one of the few lecithotrophically developing nudibranchs 
found in temperate waters in the Northeast Pacific Ocean (Goddard 2004), and its recent 
range expansion may reflect long-term changes in regional nearshore circulation regimes 
and productivity (Rebstock 2003). These factors make P. hiltoni a particularly interesting 
system for studying the complex interactions between climate change, range shifts, and 
marine invasion biology. 

If the dispersal of P. hiltoni has been facilitated by changes of oceanographic regimes 
and warming waters at the northern edge of its range, the most likely source for the new 
populations north of Monterey Bay are Central California populations. However, it is also 
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possible that individuals from father south were introduced into northern California either 
by larval dispersal or human activities. The lecithotrophic larval development of P. hiltoni 
makes it an ideal candidate for ballast water dispersal, and less likely to disperse long 
distances naturally in response to environmental changes. Two of the busiest commercial 
ports in North America are located in California, with the Los Angeles-Long Beach Har¬ 
bor well within the historic range of P. hiltoni , and the Oakland-San Francisco Harbor 
(San Francisco Bay) in the center of the extended portion of the range. However, the ab¬ 
sence of P. hiltoni from San Francisco Bay along with its prevalence in open-coast rocky 
reefs (Goddard et al. 2011) contradicts the ballast water introduction hypothesis. Another 
possible vector for the spread of P hiltoni could be small vessel traffic between regional 
ports and harbors (Wasson et al. 2001), but the mobile hunting behavior of this species 
makes it an unlikely fouling organism. Regardless of the mechanism of dispersal, Schultz 
et al. (2011) noted that R hiltoni has persisted at higher latitudes despite ocean temperature 
fluctuations from El Nino/La Nina cycles, suggesting this species may be an indicator of 
faunal range shifts due to climate change. 

Although the range expansion and some of the associated ecological effects of R hiltoni 
are well documented (Goddard et al. 2011), many questions remain. No genetic studies 
have been conducted on P. hiltoni , thus the population structure of the species is unknown, 
hampering our ability to understand the mechanisms of dispersal. Also, very little atten¬ 
tion has been paid to the southern range limit of P. hiltoni , leaving unanswered the question 
as to whether recent observations indicate a northern range expansion or overall range 
shift. In the present study, we examine population structure in R hiltoni , explore the ge¬ 
netic signature of its range expansion, and hypothesize possible dispersal mechanisms into 
Northern California. Additionally, we reviewed historical collection data near the south¬ 
ern range of P. hiltoni in order to document population density changes that may help to 
understand the population dynamics of this species. 


Materials and Methods 

Whole specimens and tissue samples of Phidiana hiltoni (Table SI) were obtained from 
different sources and various locations along the California coast (Fig. 1). Some individu¬ 
als were collected at the shoreline during low tide, other specimens and tissue samples were 
provided by colleagues or obtained from collections of the Natural History Museum of Los 
Angeles County (LACM), the Santa Barbara Museum of Natural History (SBMNH) and 
the California Academy of Sciences (CASIZ). Fieldwork was conducted under the Cal¬ 
ifornia Department of Fish and Wildlife permit #13256. Specimens collected in the field 
were preserved in 95% ethanol and deposited at the California State Polytechnic University 
Invertebrate Collection (CPIC). 

DNA was extracted from sixty specimens (Table 1) using a DNeasy Blood and Tis¬ 
sue Kit (Qiagen, Valencia, CA) using standard protocols provided by the manufacturer. 
A fragment of the cytochrome c oxidase subunit I (COI) mitochondrial gene was am¬ 
plified and sequenced using universal primers (LCO1490 5'-GGTCAACAAATCATAAA 
GATATTGG-3', HC02198 5 / -TAAACTTCAGGGTGACCAAAAAATCA-3 / ) (Folmer 
et al. 1994). PCR reaction conditions were as follows: denaturation at 95°C for 3 min, 
followed by 35 cycles of 94°C for 45 sec, 50°C for 45 sec, and 72°C for 2 min with fi¬ 
nal elongation at 72°C for 10 min. Successful DNA amplification was confirmed using an 
agarose gel electrophoresis and ethidium bromide. The PCR products were purified with a 
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Fig. 1. Map of the historic and expanded range of Phidiana hiltoni along the western coast of North 
America. Sampling localities are indicated by arrows and coded with different grey tones by geographic 
regions. 


GeneJET PCR Purification Kit (Fermentas, Waltham, MA) using standard protocols and 
were sent for sequencing to Source Bioscience Inc. (Santa Fe Springs, CA). 

Sequences were assembled and aligned using Geneious v8.1.8 (Kearse et al. 2012). The 
geographic distribution of mtDNA haplotypes was visualized by producing a haplotype 
network using the program PopArt vl.7 (Leigh and Bryant 2015) using the TCS option. 
Haplotypes were pattern-coded by locality. Genetic structure within and among popula¬ 
tions and among groups was examined using analysis of molecular variance (AMOVA) as 
implemented in Arlequin v3.5 (Excoffier and Lischer 2010). Three different AMOVA anal¬ 
yses were run to test for the effects of arranging populations into different groups based 
on the distribution of genetic variation. In the first AMOVA test, populations established 
after 1977 were included in the Northern California group while historic range populations 
were divided into two groups: Central California (populations north of Point Conception) 
and Southern California (populations south of Point Conception); this is the hypothesized 
biogeographic structure if Point Conception acts as a biogeographic barrier (Blanchette 
et al. 2008). To examine possible genetic similarities between newly formed populations 
in Northern California and those from southern Monterey Bay, two additional AMOVAs 
were run with different group arrangements, to examine whether this resulted in different 
distributions of genetic diversity among groups and among populations within groups. In 
the second AMOVA, the groups were kept the same except central California populations 
were split into two groups, Northern Central California (southern Monterey Bay) and 
Southern Central California (populations further south), see Blanchette et al. (2008). In 
the third AMOVA populations from southern Monterey Bay were pooled with Northern 
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California. Significance of the AMOVAs was tested using 16,000 permutations of individ¬ 
uals between groups. Arlequin v3.5 was also used to calculate pairwise <J>st between pop¬ 
ulations (1,000 permutations). Because the Southern California populations (Naples and 
San Clemente) were represented by one sequence each, populations were pooled together. 

Microsatellite loci were identified through sequencing by synthesis with a MiSeq plat¬ 
form (Illumina, Inc., San Diego, CA). DNA was extracted from a single specimen collected 
from Pillar Point, California (CASIZ 190249), tagged with a unique barcode during library 
preparation, and pooled with other samples for Illumina sequencing. Sequencing was con¬ 
ducted at the UCLA Genotyping and Sequencing Core facility. Automated screening of 
sequences for tetranucleotide repeats and primer design were performed simultaneously in 
MSATCOMMANDER vl.0.8 (Faircloth 2008). Twenty-two primer pairs were purchased 
from Eurofins (Louisville, KY) with a M13 tail added to the 3' end of each forward primer 
sequence. 

Five of the twenty-two primer pairs were tested with ten specimens that consistently am¬ 
plified for mtDNA to determine the PCR protocol. PCR protocol settings for the primers 
were optimized from a standard protocol by adjusting the annealing temperatures and 
elongation times until amplification was achieved. The PCR Master Mix for each locus in 
these tests included the forward primer with a M13 tail, reverse primer, and BSA (bovine 
serum albumin) and used Thermo Fisher Platinum Hot Start PCR Master Mix. The final 
PCR reaction conditions were as follows: denaturation at 95°C for 15 min, followed by 
30 cycles of 94°C for 30 sec, 60-65°C for 30 sec, and 68°C for 30 sec with final elonga¬ 
tion at 60°C for 10 min. All twenty-two primer pairs were tested with ten specimens that 
consistently amplified for mtDNA using the above conditions. Of the twenty-two primer 
pairs tested, ten polymorphic loci amplified reliably. Using the ten reliable primer pairs 
(Table S2) and the above amplification conditions, PCR was carried out with fifty-seven 
specimens. The PCR Master Mix for each locus in this final round now included a fluores¬ 
cent Ml3 tag (5'-[6-FAM] AGGGTTTTCCCAGTCACGACGTT-3') along with the orig¬ 
inal components. Genotyping was outsourced to Laragen Incorporated (Culver City, CA). 
Genotypes were scored using the Microsatellite Analysis External Plugin vl.4.4 imple¬ 
mented in Geneious v8.1.8 using the Two Surrounding Peaks setting (Kearse et al. 2012). 

In total, fifty-five individual specimens were genotyped for all 10 microsatellite loci. This 
is a small sample size for this type of study, but specimens were difficult to obtain in the 
field. Collecting sea slugs is serendipitous in nature and after two years of fieldwork only a 
small number of specimens was obtained. Additionally, most museum specimens examined 
were unsuitable for molecular work. Population subdivision in the nuclear genome was in¬ 
ferred using STRUCTURE v2.3.4 (Pritchard et al. 2000) with the default parameters; 5 
replicates for each value of K were run for 1,000,000 MCMC iterations following a burn-in 
period of 100,000. To detect the true number of clusters ( K) using the Evanno Method 
(Evanno et al. 2005) the result file from STRUCTURE was processed with STRUCTURE 
Harvester vO.6.9.84 (Earl and vonHoldt 2012). Using the selected K value (3) the result¬ 
ing files were processed with CLUMPP vl.1.2 (Rosenberg et al. 2002) and Distruct vl.l 
(Rosenberg, 2004) to generate a graphic display of the population structure. AMOVA and 
Fst pairwise genetic differentiation comparisons between populations were conducted fol¬ 
lowing the same methodology as in the mtDNA analyses. Microsatellite data were also 
analyzed via Discriminant Analysis of Principal Components (DAPC) using the adegent 
package in R (Jombart et al. 2010). 

To determine whether the abundance of P. hiltoni in Southern California has changed 
since the mid-20 th century, counts of nudibranchs by James R. Lance dating from 1953 
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to 2001 1 at six rocky intertidal sites in San Diego County (Point Loma, Hill Street, 
False Point, Bird Rock, Windansea, and South Casa Reef) were examined and ana¬ 
lyzed. The data for P. hiltoni were extracted and grouped by site, decade, and before 
and after 1963, the year when, excepting one brief trip to Bird Rock in 1956, Lance 
started sampling outer coast sites in San Diego County other than Point Loma. Counts 
made on consecutive or near-consecutive dates at any given site were excluded from 
analysis in order to reduce autocorrelation in the data; the count retained was the one 
with the highest number of P. hiltoni. Twenty-one additional counts conducted by ei¬ 
ther JG or CK from 2000 to 2016 at 4 of the same sites (Point Loma, Hill Street, 
Bird Rock, and South Casa Reef) were also included in the analysis. A Wilcoxon 
sign-rank test was implemented in JMP vl3, SAS Institute Inc. (Cary, NC) and used 
to compare the number of P. hiltoni found at Point Loma before and after 1963. 
Additional information on the recent occurrence of P. hiltoni in San Diego County 
was obtained from the website iNaturalist (https://www.inaturalist.org/taxa/48724- 
Phidiana-hiltoni) and the species database on Divebums, a San Diego dive website 
(http://species.divebums.com/index.php?l=sciname&n=Phidiana%20hiltoni), and con¬ 
firmed by the authors. 


Results 

The haplotype network of the mitochondrial COI gene recovered five distinct haplo- 
types (Fig. 2). Thirty-three individuals spanning all seven populations share the most com¬ 
mon haplotype. Three haplotypes diverge from the most common haplotype by only two 
nucleotides. The most common of these three haplotypes was found exclusively in four¬ 
teen individuals from Pillar Point, Northern California. The other two haplotypes are only 
found in specimens originating from Cayucos and Cave Landing, Central California. An 
additional specimen from Cave Landing possessed a haplotype diverging from the most 
common haplotype by three nucleotides. 

In the first AMOVA test (Northern California: Pillar Point; Central California: Cave 
Landing, Jalama Beach, Cayucos, and Carmel Point) most of the genetic variation is re¬ 
covered among groups (70.05%) and within populations (30.59%), with virtually no vari¬ 
ation among populations within groups (-0.63%) (Table 2). In the second AMOVA test 
(Northern California: Pillar Point; Northern Central California: Carmel Point; Southern 
Central California: Jalama Beach, Cave Landing, Cayucos) most of the genetic variation 
is again among groups (65.94%) and within populations (32.82%) and very little variation 
among populations within groups (1.24%) (Table 2). In the third AMOVA test (Northern 
California: Pillar Point, Carmel Point; Central California: Jalama Beach, Cave Landing, 
Cayucos) most of the genetic variation is found again among groups ( 50.47%) and within 
populations (33.5%) however genetic variation is found among populations within groups 
(16.04%) (Table 2). 

A pairwise <J>st test was run on all populations and resulted in relatively high values be¬ 
tween Pillar Point and all other populations and also between Carmel Point and Cayucos. 
The only significant difference in genetic variation found was between Pillar Point, North¬ 
ern California and each of the four Central California populations: Jalama Beach, Cave 


1 Goddard. J.H.R. 2013. Opisthobranch gastropods observed on the outer coast of San Diego County, 
California by James R. Lance, 1953-2001. knb.298.2. [online] California Academy of Sciences. Available 
https://knb.ecoinformatics.Org/knb/metacat/knb.298.2/knb [2017 Jun 15]. 



SOUTHERN CALIFORNIA ACADEMY OF SCIENCES 


.5 H 

*oT s 

<=> £ 
o J3 

vi u 

Ci 3 
^ 3 

^ 00 
<L> ^ 

_£ 3 


3 05 

o r 2 
'-C U 


3 Sfe 

>, e 

U 1 

ox J 

<L> 

I 8 

j u 

<D -s 
£ % 
U & 

•S g 

3 S 

<D 3 

OQ ^ 


.2? 53 

C/5 

r 4 —* 

*05 3 

rd O 

> 00 

c ^ 

• — g 

-* 3 

a" c 

jl> C" 

Vh r-< 

^ o 

j 5 Cl 

’> E 

^ s 

’"O *-< 

<D CO 

.£ U 

5 

X) o 
O w 

V5 ^ 

w> 3 
.2 U 

CL 

3 OX 

O G 4-^ G 

ox *"3 .5 "o 
^ c o Cl 
Cl . 


3 CO 

E ’£ 

3 ut 

o,o 

•g 1 

£ u 

«§ 13 

OO S-N 

co 'S 

u g 

3 U 

S—I *■ 


<u .3 

U 

c &- 

<3 13 




° U 


.2 J 

_3 ^ 

3 CO 

1° 

Qh 

<u O 
4—> CO 

co <u 

b CQ 

CO 

CL CO 

8 £ 

4> -2 

<D 3 


0 .2 

o c 

3 O 
3 £ 
"O 3 

D o 

CL, _ 
>> 3 

o ~ 

cl, 53 

2 U 


<L> ™ 

= £ 

3 

O .2 

,N —^ r-* 

1 I 

C ;_g 

J§ u 
u £ 

o 

co Ho 


c 

0/ 

u 


a-) 

JC 

4—► 

O 

2 


.2 

*3 

"a 

o 


£ 


c£ 

CO 


o 

2 

cd 

ox 

# c 

# S- 

3 

2 ~ 
o E 

CL 

CO 

S 2 


cl 

CO 

2 


Q 


— co Cm 


o 

o 


— 

£ 

cO 


O 55 


■r o ^ L 


c 2 c 
i2 

3 3 

^ U 
o w 


O O 
<L> 

i 5 

> z 

o r. 
S r r 
< g 

'a* 

o) § 

(D 

3 


a 


^ 3 

o 

jo 


C w s 

j§ 0 £ 

n § •§ 

5 .2 ~ 

JO G co 

ego 

O 

^ *3 53 

N 0» 6 

g> £ O 

|s« 

3 7 O) 

O 3 G 

*r O cO 

O W ^ 
-a ^ 

V5 c O 
^ cO y 

O? O 
JL CJ \J 


CD 

On 

uo 

r- 

lO) 

o\ 

3 


r- 

Gn 

o 

sO 


o 


04 

04 

On 

3 


04 

O 

On 

TM« 

> 

o 

o 

Ol 

O 

o 


o 

o o 

o 

O O 

O 


V5 

1) 

O 


"O 

Ol 

ro 


04 

Os 

DO 

^t* 

Ol 

^o 


* 

»—-* 

^J" 

"0- 


ro 

r- 

o 

^o 

* ^ 


^t 

o 

^O 


On 

ro 

105 

^J" 


04 

ON 

o 

ro 

o- 

105 

Ol 

^O 

o 


o 

VO 

r- 

O 

^o 

^O 

ro 

^o 

105 

CO 

o 

1 

o 

o 

o 

o 

O 

O 

o 

o 


3 

E 


co 

> 

t4-i 

o 

Ln 


a 

o 

CL 

£ 

O 

o 

(U 

o 

G 

.2 


V5 

(U 


CT 

C-^H 

O 

£ 

23 

00 


.o 

4—> 

.2 

CO 

> 


CD 

CJ 


o 

00 


r] \j- ^ i^ 
"0" ^ 


a 

E 

4—> 


n n ^ 

"0" "0" 


O 

E 

*1 


cn 


o o 


•o 

m 

On. 


"0" 

04 


"0” 


o 

VO 

»o 

Ov 

Oj 

oo 


o 

105 

o 

O 

o 

»o 

—.’ 

oi 

o 

vd 

CO 


i 

ro 

VO 


co 

105 


CO 


Gn 


*05 

— 

Ol 

*05 

*05 

04 

co 


*05 

oo 

Ol 

o 

OO 


o 

*05 

oc 

Of 

Ol 

CO 

OO 

CO 

VO 

co 

\o 

o 

*05 

*05 

VO 



o 

^o 

oo 

CO 

o 


oc 

GN 

o 

of 

of 

Ol 


of 

CO 

CO 

o 

— 

of; 

Ol 

o 

— 

of 

04 

o 

’—' 

Of 

o 

o 

o 

o 

O 

o 

O 

O 

o 

o 

o 

o 


oo 

CO 

— 

04 

oo 

of 


04 

o 

04 


04 

*05 

CO 

04 

—4 

o 

04 

04 


o 

Gv 

04 

—« 

04 

*05 

o 

OO 

04 

*05 

o 

OO 

Ov 

OO 

o 

OO 

OO 

o 

vd 

of 

OO 

o 

vd 

of 

vd 

—* 

vd 

of 


04 


"0- — r- 
"0" "0" 


V5 

zn 

zn 

CL 

Dh 

CL 

3 

3 

3 

O 

O 

O 


O 



zn 

3 

75 




V5 

3 

V5 

r - 1 




zn 

V5 



5 





o 




o 

^5 

VI 

CL 

E 

4^ 

7Z 



CL 

E 

3 



75 

CL 

*4—* 

E 

*4-> 

3 

3 

3 




3 

3 




3 

3 


O 

a 

CL 

c£ 

CL 

O 

CL 


04 

OX 

O 

Vh 

a 

CL 

£ 

CL 

O 

Cl 


ro 

OX 

p 

a 

CL 

t£ 

Oh 

c£ 

to 

O/j 



E 

to 

o/j 

r - 1 


( C 

ox 

ox 

f4 

r~] 

3 

2 h 


‘5L 





*S- 



i —1 

o 

O 


3 

3 

o 

o 

r-j 

3 

3 

o 

o 


g 

£ 

.G 

_G 

£ 

g 

£ 

•CC 

4—* 

_o 

o 

S-N 

£ 

£ 

4—» 


ox 

G 


3 

2 

a < < 


t2 6 < < ^ e2 


cO 

I 






RANGE SHIFT OF PHI DIANA H1LTONI 


9 


CK130 

CK131 


CK060 

CK061 

CK062 

CK063 

CK064 

CK066 

CK067 


CK068 

CK069 

CK070 

CK073 

CK074 

CK075 

CK076 


CK029 

CK030 

CK031 

CK032 

CK033 

CK077 

CK079 

CK081 


CK078 



I 10 samples } 

VQ% 

1 sample 

Northern California 

Pillar Point 

Central California 

Cave Landing 

Cayucos 

o Carmel Point 
o Jalama Beach 

Southern California 

o Naples 

San Clemente Island 


CK088 



Fig. 2. Haplotype network of COI mitochondrial sequences generated with PopArt vl.7. Each circle 
represents a unique haplotype and its area is proportional to the number of specimens sequenced with that 
haplotype. Each pattern represents the geographic origin of the individual specimens, as indicated in the 
legend. Isolate codes are indicated next to each haplotype. 


Landing, Cayucos and Carmel Point (Ost = 0.69, p = 0.03; <J>st = 0.75, p = 0.00; <I>st = 
0.69,/? = 0.00; <t>sT = 0.67,/? = 0.04 respectively) (Table 3). This suggests genetic differen¬ 
tiation between Pillar Point and the other populations. However, this result, as well as the 
lack of significant differentiation among other pairwise comparisons, must be interpreted 
with caution due to the limited number of samples from all collection sites except for Pillar 
Point and Cave Landing. 

Analysis of microsatellite data with Structure Harvester using Evanno’s method 
(Evanno et al. 2005), a maximum value of the rate of change (A K) in the log probabil¬ 
ity of data was obtained at K = 3 (Fig. 3B). These three recovered clusters are unevenly 
distributed among geographic regions with no obvious geographic subdivision (cluster 1: 
red, cluster 2: blue, cluster 3: yellow. Fig. 3A). Moreover, all individuals exhibit a non-zero 
probability of belonging to any one of the three clusters. 

All AMOVA tests with different groupings produced very similar results; the overwhelm¬ 
ing majority of the genetic variation was recovered within populations (92.32-93.61%) and 
some among populations within groups (6.51-9.98%), with virtually no variation among 
groups (-3.59-1.17%) (Table 4). Pairwise F S t comparisons produced very low values 
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Table 3. <t>sT pairwise comparison values for mitochondrial haplotype data obtained with Arlequin v3.5 
(lower triangular) and associated p values (upper triangular). Significant values (p < 0.05) in bold. 



Pillar 

Point 

Carmel 

Point 

Cayucos 

Cave 

Landing 

Jalama 

Beach 

Naples 

San 

Clemente 

Pillar Point 

— 

0.03062 

0.00019 

0.0000 

0.03537 

0.09994 

0.99994 

Carmel Point 

0.69490 

— 

0.21358 

0.32000 

0.99994 

0.99994 

0.99994 

Cayucos 

0.75345 

0.62791 

- 

0.99994 

0.99994 

0.99994 

0.99994 

Cave Landing 

0.68730 

0.29687 

-0.05381 

- 

0.99994 

0.99994 

0.99994 

Jalama Beach 

0.66771 

0.0000 

0.0000 

-0.33043 

— 

0.99994 

0.99994 

Naples 

0.61778 

0.0000 

0.0000 

0.0000 

0.0000 


0.99994 

San Clemente 

0.61778 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

- 


across the entire range suggesting little to no genetic differentiation between populations 
(Table 5). DAPC, which attempts to group individuals using a k-means clustering algo¬ 
rithm, suggests that the entire metapopulation of P. hiltoni cannot be divided into more 
than one group based on microsatellite data (Fig. 4). 

From 1953 to 1962, Jim Lance sampled for nudibranchs on the outer coast of San 
Diego County only at Point Loma, where he found at least one Phidiana hiltoni dur¬ 
ing half of his 28 trips during that period (Fig. 5A). In 1964 he began to sample addi¬ 
tional outer coast sites. Since then, only two more P. hiltoni were found at Point Loma 
(Figure 5A & B), one by Lance in July 1968 and one by JG in June 2001. Similarly, 
P. hiltoni was observed on only about 10% of the trips to each of the other five outer 
coast sites (Fig. 5A), and was found in lower numbers per trip than had been seen in 
the earlier period at Point Loma (Fig. 5B). Significantly fewer P. hiltoni were found at 
Point Loma after 1963 than before (Wilcoxon rank sum test, P = 0.002) (Figs. 5B, C). 
Except for the 1990s, when only 10 total trips were made to two sites, P. hiltoni has 
been found in San Diego County in low numbers in each of the decades since the 1960s 
(Fig. 5C). Finally, P. hiltoni has been photographed subtidally in San Diego County 
at least 10 times since 2005 (http://species.divebums.com/index.php71 = sciname&n = 
Phidiana%20hiltoni; https://www.inaturalist.org/taxa/48724-Phidiana-hiltoni). 



Fig. 3. Genetic clustering analysis for the entire data set of 57 individuals and 10 microsatellite re¬ 
gions as estimated by STRUCTURE v2.3.4. A. Genetic clustering plot for K = 3 clusters, generated with 
CLUMPP vl.1.2. Each grey tone represents a different genetic cluster. Bar graphs show average posterior 
probability of membership (y-axis) of each individual. Populations are delimited by dark vertical lines. B. 
Graph of A K = mean (| L'(/0 1 ) / sd( L( AT)) as a function of K (potential number of genetic clusters) gen¬ 
erated by STRUCTURE Harvester vO.6.9.84. The most likely number of clusters is indicated by the modal 
value, in this case K = 3. 

























Table 4. AMOVA test results for microsatellite genotype data with three separate population groupings obtained with Arlequin v3.5, significant values (p < 0.05) in 
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Table 5. Fst pairwise comparison values for microsatellite genotype data, obtained with Arlequin v3.5 
(lower triangular) and associated p values (upper triangular), significant values (p < 0.05) in bold. The two 
southern California populations of Naples and San Clemente were combined into one due to low sampling 
numbers. 



Pillar Point 

Scott Creek 

Carmel Point 

Cave Landing 

Jalama Beach 

S. California 

Pillar Point 

— 

0.93776 

0.01968 

0.64388 

0.39049 

0.35181 

Scott Creek 

0.13279 

- 

0.99994 

0.99994 

0.99994 

0.99994 

Carmel Point 

0.14426 

0.35183 

— 

0.00306 

0.00444 

0.10229 

Cave Landing 

0.01779 

0.17272 

0.15117 

— 

0.06255 

0.05986 

Jalama Beach 

0.04219 

0.23112 

0.19336 

0.05260 

— 

0.19821 

S. California 

0.12656 

0.36283 

0.30272 

0.15087 

0.15965 

- 


Discussion 

The California coast is a prime example of a region where climate change is impact¬ 
ing native marine ecosystems through changes in ocean temperatures, seawater chemistry, 
and coastal current regimes (Barry et al. 1995; Sagarin et al. 1999; Harley et al. 2006). 
With ocean temperatures increasing, species are predicted to shift their ranges poleward, a 
trend that has already been observed across a wide range of taxonomic groups, mostly at 
temperate latitudes (Dawson et al. 2010; VanDerWal et al. 2012). It is difficult to predict 
how environmental change, including biotic exchanges resulting from species range shifts 
and introductions, will affect ecological systems. Heilman et al. (2008) emphasized that 



Fig. 4. Results of the DAPC analysis indicating that the data, when analyzed as principal components, 
cannot be divided into more than one group. 
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Fig. 5. Occurrence of Pliidiciua hiltoui at six rocky intertidal sites in San Diego County, 1953-2016. A. 
Frequency of occurrence by site, 1953-1962 and 1964-2016. Number of trips to each site shown above bars. 
B. Number of P. liiltoni found on each trip, 1953-2016. C. Mean number (± SE) of P. hiltoui found per trip 
at Point Loma and other sites combined, by decade. 



















14 


SOUTHERN CALIFORNIA ACADEMY OF SCIENCES 


global climate change compounds this difficulty because it transforms transport and in¬ 
troduction mechanisms, impacts the distribution of existing invasive species, and alters the 
effectiveness of control strategies. Efforts to manage and conserve marine ecosystems in the 
face of climate change will require improvements to the existing predictive framework to 
aid in preventing future introductions (Harley et al. 2006). In this context, understanding 
the mechanisms behind range expansions (or shifts) of individual species will contribute 
to a larger body of evidence, critically important for predicting the biological effects of 
climate change. 

The data presented in this paper provide insight into the processes underlying the range 
expansion in P. hiltoni. The analysis of mitochondrial DNA (mtDNA) sequence data of 
Pliidicma hiltoni across both its historical range and extended range in Northern California 
revealed limited genetic structure. The COI haplotype network shows very little polymor¬ 
phism with only five haplotypes in total being recovered (Fig. 2). The most common of 
these haplotypes is shared among individuals from all seven populations, including the 
newly formed populations in Northern California. However, there are a few haplotypes 
only detected in certain populations. For example, one of the two haplotypes found at Pil¬ 
lar Point, north of Monterey Bay, was found nowhere else. Because populations north of 
Monterey Bay did not exist prior to 1977 (Goddard et al. 2011), it is likely that this unique 
haplotype is also present south of Monterey Bay but at such low frequency that it has not 
been yet detected, and may have become more common in the recently colonized popu¬ 
lations due to founder effects. To better understand the geographic structure of P. hiltoni 
based on mtDNA, three AMOVA analyses were run, each with different population group¬ 
ings. In each of the three groupings, the highest percent variation was consistently found 
among groups (Table 2). However, there is a decrease in this percent variation as the group¬ 
ings structure are altered (70.05%, 65.94%, 50.47% respectively), suggesting that the first 
grouping arrangement (in which the Northern California group includes only populations 
found in the extended range) best represents population genetic structure according to 
mtDNA data. Pairwise <4>st comparisons agree with the geographic structure recovered in 
the AMOVAs (Table 3). Relatively high <J>st values indicating genetic differentiation were 
found between Pillar Point and the Central and Southern California populations, and be¬ 
tween Carmel Point and the two populations of Cayucos and Cave Landing. This is con¬ 
sistent with the haplotype network results, showing that Pillar Point, Cayucos and Cave 
Landing possess divergent haplotypes from the most common haplotype found across the 
range of P hiltoni . 

While analyses of mtDNA sequences suggest population structure and genetic differen¬ 
tiation among groups, no genetic pattern corresponding to geography was detected using 
microsatellite data. Across the range of P. hiltoni , Structure analyses showed several indi¬ 
viduals have nearly identical probabilities of cluster membership, a pattern that is partic¬ 
ularly apparent at Cave Landing. This genetic uniformity could be an indication of self¬ 
recruitment at Cave Landing, where the concavity in the coastline (accentuated by the 
750-meter long rock jetty on Point San Luis) may encourage larval retention. Goddard 
et al. (2011) suggested that population structure in P hiltoni should be affected by up- 
welling shadows, resulting in local retention of short-lived larvae, particularly at the north¬ 
ern end of bights along the coastline (Graham and Largier 1997; Roughan et al. 2005). 
AMOVAs on microsatellite data included the same three distinct groupings as mitochon¬ 
drial analyses except for the addition of Scott Creek to the Northern California group and 
the removal of Cayucos (where microsatellite data were not successfully recovered). In all 
AMOVA tests, most genetic variance was found within populations (93.61%, 92.32% and 
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93. 11%, respectively) (Table 4). Consistently low variation among groups is indicative of 
high levels of gene flow across the range, as seen in the mtDNA data. Pairwise compar¬ 
isons of Fst values from microsatellite data are also consistent with high gene flow across 
the range, with generally low Fst values between populations (Table 5). The slightly higher 
Fst values found between Scott Creek and Carmel Point, as well as between Carmel Point 
and Southern California are unreliable due to the low sample sizes from Scott Creek and 
Southern California. 

AMOVA and pairwise $st comparisons using mtDNA data suggest Northern Califor¬ 
nia populations are genetically distinct, and consequently the origin of the recently founded 
populations remains unclear. On the contrary, microsatellite data indicate very high levels 
of gene flow in P. hiltoni, with Northern California populations genetically indistinguish¬ 
able from those in Central California. This discrepancy of results from mtDNA vs. nuclear 
data might be explained by the reproductive behavior of P. hi/toni. Rutowksi (1983) found 
that species with high rates of cannibalism require several couplings, usually with differ¬ 
ent mates, in order for all the eggs in the egg mass to be fertilized. It is very likely that 
R hiltoni (considering the pugnacious and cannibalistic nature of this species) exhibits ab¬ 
breviated coupling times, requiring several mates to fully fertilize egg masses. Additionally, 
P. hiltoni exhibits locally high population densities in Central California, potentially facili¬ 
tating promiscuity and thus intense sperm competition. The main consequence of this mat¬ 
ing system would be that egg masses produced by a single individual (functional female) 
will be sired by multiple partners, generating offspring with genetically identical mitochon¬ 
dria but different nuclear alleles. If recently established populations in Northern California 
are the result of sporadic events involving dispersal of larvae (or rafting of egg masses) 
produced by a small number of females, it is likely that founder effects facilitated retention 
of these rare mitochondrial haplotypes. While the diversity of nuclear alleles should also 
decrease due to drift and founder effects, the polyandrous mating system of R hiltoni has 
the potential to mitigate these effects. Under this scenario, Northern California popula¬ 
tions could harbor mitochondrial haplotypes that were previously very rare in the historic 
range, while nuclear alleles from the source population are more broadly represented. If 
this hypothesis is correct, additional sampling from Central California should detect all 
or most Northern California haplotypes. An alternative explanation is that the lack of ge¬ 
netic structure in microsatellite data is an artifact of the limited sample size. If this is the 
case, additional sampling across the range of R hiltoni would improve the reliability of the 
results of this study. 

Another outstanding question is what mechanism(s) allowed R hiltoni to cross Monterey 
Bay starting in the late 1970s, or what physical or biological barriers restricted the prior 
range of this species. If recently established populations in Northern California are indeed 
the result of sporadic dispersal by a limited number of individuals, this would suggest that 
this dispersal was not triggered by a gradual process, such as increasing ocean temper¬ 
atures, but instead by the temporary or intermittent opening of a corridor. One distinct 
possibility is that a weakening of the upwelling shadow in Monterey Bay (Graham and 
Largier 1997) due to climate change may have facilitated this process. Pennington et al. 
(2000) documented decadal-scale changes in the oceanographic conditions near the center 
of Monterey Bay region consistent with those described for the 1976-77 climate shift in 
the North Pacific Ocean. These include increased stratification of surface waters, warmer, 
less productive waters during non-upwclling seasons and a later onset of upwelling. These 
changes may have reduced larval retention in northern Monterey Bay, allowing R hiltoni 
to disperse into Northern California. 
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Few studies have examined the population genetics of other marine invertebrate taxa 
that have experienced recent range expansions or range shifts. In California, a similar study 
by Dawson et al. (2010) examined three hypotheses/scenarios that explain the causes of 
range limits of species, and concluded that in the volcano barnacle, Tetvaclita rubescens , the 
northern range boundary is maintained by migration load arising from flow of maladapted 
alleles into peripheral locations. Dawson et al. (2010) proposed that in this species (with 
planktonic-dispersing larvae), environmental amelioration, likely due to climate change, 
resulted in a reduction of the strength of selection against immigrant phenotypes in the 
northern range boundary, allowing the species to expand northward. The case of P. hiltoni 
is very different in several respects, but the main difference is that whereas T. rubescens was 
declining near its northern range limit, P. hiltoni has been and remains common. Framing 
of our data in the three scenarios/hypotheses proposed by Dawson et al. (2010) suggest 
physical barriers to dispersal is the most likely mechanism that historically restricted mi¬ 
gration in P. hiltoni. 

This study included a sample of 57 individual specimens collected across the range of 
P. hiltoni , but the sample size from Southern California, south of Point Conception, is 
small. Only two individuals were collected despite a substantial collecting effort by the se¬ 
nior author in this region during two consecutive years. This suggests that P. hiltoni could 
have become rare in the southern portion of its range. However, this assumption should 
be interpreted with caution. Bates et al. (2015) shown that abundance-related species de¬ 
tectability, particularly important in uncommon, difficult-to-detect marine species such as 
P hiltoni , has the potential to confound our understanding of the true location of range 
edges. Bates et al. (2015) emphasized the importance of simulation and modeling, but also 
long-term monitoring with consistent sampling effort through time. In this case, we ana¬ 
lyzed high-quality, long-term observational data collected at fairly regular intervals from 
the same region. These data suggest that the historical abundance of P. hiltoni in San Diego 
County, especially Point Loma (the type locality of P. pugnax Lance 1961 [ = P. hiltoni]) 
has declined (Figs. 4A-C). Reasons for this decline remain unknown and warrant fur¬ 
ther investigation, although rising ocean temperatures appears to be a viable hypothesis. 
Notably, the abundance of another species of sea slug, Feliniave californiensis, was once- 
common in Southern California, but became extinct there in the 1980’s (Goddard et al. 
2013). Although individuals of F. californiensis reappeared in 2003 and the species has 
since been found in a few isolated localities in Southern California (Goddard et al. 2013; 
Hoover 2015), its populations have not completely recovered. It is unclear whether there is 
a link between the decline of these two ecologically distinct species, but if there is, it may 
be a symptom of larger and more pervasive environmental change. The apparent decline 
of P hiltoni in Southern California along with its dispersal northward needs to be sub¬ 
stantiated with further monitoring and additional data analyses (Bates et al. 2015), but if 
confirmed, would suggest this is a true poleward range expansion rather than a temporary 
shift (Parmesan et al. 1999). Understanding the process by which P. hiltoni migrated north¬ 
ward may provide insight as to how other benthic organisms will respond to rising ocean 
temperatures and changes in ocean current systems (McGowan et al. 1998). 
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Appendix 


Table SI. Complete list of specimens sequenced for this study, including isolate number, locality, collec¬ 
tion date, and GenBank accession numbers. 


Isolate 

Locality 

Collection date 

GenBank accession # 

CK11 

San Clemente Island, CA 

1961 

MK333330 

CK29 

Cave Landing, CA 

10/8/2014 

MK333291 

CK30 

Cave Landing, CA 

10/8/2014 

MK333292 

CK31 

Cave Landing, CA 

10/8/2014 

MK333293 

CK32 

Cave Landing, CA 

10/8/2014 

MK333294 

CK33 

Cave Landing, CA 

10/8/2014 

MK333295 

CK34 

Naples, CA 

12/2009 

MK333313 

CK46 

Tarantula Reef, Jalama Beach, CA 

12/15/2009 

MK333331 

CK47 

Tarantula Reef, Jalama Beach, CA 

12/15/2009 

MK333332 

CK48 

Tarantula Reef, Jalama Beach, CA 

12/15/2009 

MK333333 

CK49 

Tarantula Reef, Jalama Beach, CA 

12/15/2009 

MK333334 

CK50 

Tarantula Reef, Jalama Beach, CA 

12/15/2009 

MK333335 

CK51 

Tarantula Reef, Jalama Beach, CA 

12/15/2009 

MK333336 

CK52 

Tarantula Reef, Jalama Beach, CA 

12/15/2009 

MK333337 

CK53 

Tarantula Reef, Jalama Beach, CA 

12/15/2009 

MK333338 

CK60 

Pillar Point, CA 

6/18/2015 

MK333314 

CK61 

Pillar Point, CA 

6/18/2015 

MK333315 

CK62 

Pillar Point, CA 

6/18/2015 

MK333316 

CK63 

Pillar Point, CA 

6/18/2015 

MK333317 

CK64 

Pillar Point, CA 

6/22/2015 

MK333318 

CK65 

Pillar Point, CA 

6/22/2015 

MK333319 

CK66 

Pillar Point, CA 

6/22/2015 

MK333320 

CK67 

Pillar Point, CA 

6/22/2015 

MK333321 

CK68 

Pillar Point, CA 

6/22/2015 

MK333322 

CK69 

Pillar Point, CA 

6/22/2015 

MK333323 

CK70 

Pillar Point, CA 

6/22/2015 

MK333324 

CK71 

Pillar Point, CA 

6/22/2015 

MK333325 

CK73 

Pillar Point, CA 

6/22/2015 

MK333326 

CK74 

Pillar Point, CA 

6/22/2015 

MK333327 

CK75 

Pillar Point, CA 

6/22/2015 

MK333328 

CK76 

Pillar Point, CA 

6/22/2015 

MK333329 

CK77 

Cave Landing, CA 

1/16/2016 

MK333296 

CK78 

Cave Landing, CA 

1/16/2016 

MK333297 

CK79 

Cave Landing, CA 

1/16/2016 

MK333298 

CK80 

Cave Landing, CA 

1/16/2016 

MK333299 

CK81 

Cave Landing, CA 

1/16/2016 

MK333300 

CK82 

Cave Landing, CA 

1/16/2016 

MK333301 

CK83 

Cave Landing, CA 

1/16/2016 

MK333302 

CK84 

Cave Landing, CA 

1/16/2016 

MK333303 

CK86 

Cave Landing, CA 

1/16/2016 

MK333304 

CK87 

Cave Landing, CA 

1/16/2016 

MK333305 

CK88 

Cave Landing, CA 

1/16/2016 

MK333306 

CK89 

Cave Landing, CA 

2/5/2016 

MK333307 

CK93 

Cave Landing, CA 

2/5/2016 

MK333308 

CK94 

Cave Landing, CA 

2/5/2016 

MK333309 

CK98 

Cave Landing, CA 

2/5/2016 

MK333310 

CK117 

Carmel Pt, Monterey Bay CA 

— 

MK333289 

CK118 

Carmel Pt, Monterey Bay CA 

— 

MK333290 

CK130 

Cayucos, CA 

— 

MK333311 

CK131 

Cayucos, CA 

— 

MK333312 
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Table S2. List of primer pairs (and their sequences) used to amplify polymorphic microsatellite loci in 
P. hiltoni. Bolded portion of forward primers indicate M13 tail. 


Primer 

Sequence 

Phil760625F 
Phil760625R 
Phil792112F 
Phil792112R 
Phil820905F 
Phil820905R 
Phil928092F 
Phil928092R 
Phil98151F 
Phil98151R 
Phill09255F 
Phill09255R 
Phil585958F 
Phil585958R 
Phil918696F 
Phil918696R 
Phil 121774F 
Phil 121774R 
Phil315595F 
Phil315595R 

AGGGTTTTCCCAGTCACGACGTTAACGTCGTCATGGAATTCACAG 

GTTTATTAATGGCGGCGATGTGAC 

AGGGTTTTCCCAGTCACGACGTTAACCAATCGACGACAAGCTAAC 

GTTTGTCTCCGTGTTAAGTGTTGC 

AGGGTTTTCCC AGTC ACG ACGTT AC ATTACT CC ACT CG ACT CAGG 
GTTTAGTCTCGGTCCATGAATCAGG 

AGGGTTTTCCC AGTCACGACGTTGATTCTATGCCACACACCTTGG 
GTTTAATGTATCTGCTTCATCCGTGC 

AGGGTTTTCCCAGTCACGACGTTAGAGGAATAGTCGCGGAACTAC 

GTTTCATCATTGCGTCAGATGTCC 

AGGGTTTTCCCAGTCACGACGTTCACACGTTCATACACTCACCTG 

GTTTAACACCGAGACAAGACATGC 

AGGGTTTTCCC AGTCACGACGTTACTCTCTCACACCTGTCAAGTC 
GTTTCACCTCAGTACAGTCTCGTG 

AGGGTTTTCCC AGTC ACG ACGTTACTCTCTCACACCTGTCAAGTC 
GTTTCACCTCAGTACAGTCTCGTG 

AGGGTTTTCCC AGTC ACG ACGTTGTC A AGTG A ATA AG ACGGCG AG 
GTTT CTGCCT GCTATAC AT CC ATCC 

AGGGTTTTCCC AGTCACGACGTTGTAAC AC AGTGTCCGTATGTGG 
GTTTATCATTCTACGTGCATGCTGTC 






